Compute stacking weights with an active-set Newton method - #412

Fixes #119.

This PR was assisted by Claude.

Summary

stacking_weights() optimized the weights with stats::constrOptim(), a log-barrier method with an inner BFGS. It is slow for many models, as #119 reports (more than 10 minutes for 10 000 models). It also often stops far from the optimum: constrOptim() reports convergence after 2 or 3 outer iterations, and the barrier term keeps the weights away from zero. loo gave no warning.

This PR finds the weights with an active-set Newton method that uses only base R and matrixStats.

@yao-yl @jgabry, the stacking weights change for many inputs, so I would like your review.

How far master was from the optimum

For stacking with fixed predictive densities, a bound shows how far any weights are from the optimum. With p_i = sum_k w_k exp(lpd_ik) and g_k = mean_i exp(lpd_ik) / p_i, Jensen's inequality gives that the log score of w is at most N * log(max_k g_k) below its maximum. At the optimum, g_k = 1 for every model with positive weight and g_k <= 1 for the others (Lindsay, 1983). The same bound appears for log-optimal portfolios (Cover, 1984), which is the same optimization problem.

I simulated pointwise log predictive densities for N observations and K normal models, 50 data sets per cell, and computed this bound for the master weights:

The gap is in units of the log score. A weight difference near 1 means that master put the weight on a different model than the optimum. Larger iteration limits for constrOptim() reduce the gap but do not remove it.

The old code could also give a non-zero weight to a clearly worse model. The optimization starts from equal weights, and when it stops early, models whose optimal weight is zero keep part of their starting weight. With 10 models, master gave more than 0.01 weight to a model with optimal weight zero in 54 of 100 simulated data sets for N = 100, and in 70 of 100 for N = 1000. In 27 and 57 of these, that model's elpd was more than 2 SE below the best model, and the total weight on such models was up to 0.89. With 3 or 5 models this happened in 2 of 400 data sets.

Method

Stacking with fixed predictive densities is the maximum likelihood estimate of mixture proportions with fixed components. stacking_weights() now calls the internal stacking_weights_active_set(), a variant of the constrained Newton method (CNM) of Wang (2007):

1. Start with the model that has the largest total log predictive density, and add the models needed so that every observation has a finite predictive density.
2. Optimize the weights of the active models with Newton steps on the simplex (a KKT system for the sum constraint, with a line search). A model whose weight reaches zero leaves the set.
3. Compute g_k for all models. Stop when max_k log(g_k) is at most 1e-8, so the log score is at most N * 1e-8 below its maximum. Otherwise add the model with the largest g_k and go to step 2.

CNM adds all local maxima of the gradient function at each iteration and solves the quadratic approximation with a least-squares formulation. This implementation adds one model per iteration and solves the KKT system directly, which is enough here because the optimal weights are usually sparse. Each check over all models costs one pass over the N x K matrix.

Over 1800 simulated data sets (N from 20 to 1000, K from 2 to 50), the largest gap bound with this PR is 9.3e-7. Two other approaches were not good enough. A softmax parametrization with optim(method = "L-BFGS-B") is fast but stalls with gaps from 0.1 to 3, because the gradient of a weight is proportional to the weight. EM for mixture weights (Cover's multiplicative update), also with SQUAREM acceleration, reaches the optimum but took more than 60 s for 1000 models.

?loo_model_weights describes the method in Details and cites Wang (2007) and Lindsay (1983).

References:

• Wang, Y. (2007). On fast computation of the non-parametric maximum likelihood estimate of a mixing distribution. Journal of the Royal Statistical Society: Series B, 69(2), 185–198.
• Lindsay, B. G. (1983). The geometry of mixture likelihoods: a general theory. The Annals of Statistics, 11(1), 86–94.
• Cover, T. M. (1984). An algorithm for maximizing expected log investment return. IEEE Transactions on Information Theory, 30(2), 369–373.

API change

optim_method and optim_control of stacking_weights() and loo_model_weights() no longer have an effect. Their defaults are now NULL. If a user supplies them, a warning says that they are deprecated and ignored. They can be removed in a later release.

The stacking weights change for many inputs, also with 3 models, because the new weights are at the optimum. For the 3-model test data in test_model_weighting.R, the weights change from 0.203658, 0.796342, 3.6e-8 (gap bound 3.8e-4) to 0.203608, 0.796392, 0 (gap bound 2.4e-13). Please add a NEWS entry; I suggest a text below.

exp_diff_over_exp() in R/helpers.R computed the gradient of the old objective and has no caller after this change. This PR removes it and its tests.

Tests

tests/testthat/test_model_weighting.R:

• stacking_weights gives expected result compared with values from constrOptim(). It now checks the gap bound and compares with an EM solution of the same data.
• New: 500 models, gap bound below N * 1e-8 and sparse weights.
• New: observations that only a combination of three models can predict (-Inf entries), a duplicated model, and identical models.
• New: the deprecation warning in stacking_weights() and loo_model_weights(), and the same weights with and without the deprecated arguments.
• The serialized snapshot of the 3-model weights has the new values.

With the master R/loo_model_weights.R, the new and changed tests give 10 failed expectations. With this PR, devtools::test() with NOT_CRAN=true passes all 1269 tests and gives no warnings (R 4.6.1, testthat 3.3.2, matrixStats 1.5.0, Linux x86_64).